This document is a follow up to one prepared for an initial briefing to Dr. Angell on analyzing the equity implications of Public Safety Power Shutoffs (PSPS) in California.

This analysis mirrors the earlier one except that is uses the actual areas impacted by PG&E PSPS in 2019, rather than Fire Hazard Severity Zones, with two (2) statewide datasets describing health equity:

Data on the 2019 PG&E Public Safety Power Shutoff Events

Spatial data of PSPS affected areas initially received from PG&E was cross-referenced with reports on their website.

There are some reports on PG&E’s website for which we do not have spatial data about affected areas.

Spatial data from PG&E contained boundaries for 12 separate areas, which - using layer names - were cross-referenced to reports for 7 specific de-energization events (shown on the map below, click for PG&E’s estimates start/end times and number of affected customers).

event1 <- filter(pge_events, event == "June_7-9")
event2 <- filter(pge_events, event == "September_25-27")
event3 <- filter(pge_events, event == "October_5")
event4 <- filter(pge_events, event == "October_9-12")
event5 <- filter(pge_events, event == "October_23-25")
event6 <- filter(pge_events, event == "October_26-November_1")
event7 <- filter(pge_events, event == "November_20-21")

pal2 <- colorQuantile(palette = c("#2172B4","#9DCAE1","#A4DBA1","#008836"), n=4,  domain = 0:100)


maphpi <- merge(tracts_sp, { 
          
          hpi[,.(ct10, hpi2_pctile)]
  
          
        })

    leaflet() %>%
     # addProviderTiles(providers$Esri.WorldStreetMap, group = "Street Map") %>%
     #  addProviderTiles(providers$Stamen.Terrain, group = "Terrain") %>%
     #  addProviderTiles(providers$Esri.WorldImagery, group = "Sattelite") %>%
    addProviderTiles(providers$Stamen.Toner, group = "B/W") %>%  
        addPolygons(data = maphpi,
          
        color = "#444444",
        weight = 1,
        smoothFactor = 0.1,
        fillOpacity = 0.5,
        fillColor = ~ pal2(maphpi$hpi2_pctile),
        stroke = FALSE,
        popup = paste0("This location is in the top ",round(maphpi$hpi2_pctile,0),"the percentile of all tracts in the state for the healthy places index"),
        group = "HPI"
      ) %>%
    addPolygons(
      data = st_zm(event1),
      color = "#252525",
      weight = 2,
      smoothFactor = 0.1,
      fillOpacity = 0.7,
      fillColor = "#377eb8",
      stroke = TRUE,
      label = ~event1$event,
      popup = paste0("This PSPS event spanned ", round(event1$Duration_hrs,0), " hours   \n   (beginning ", as.character(event1$start_time), " and ending at ", event1$end_time,").     \n   Ultimately PG&E estimates that it affected ",format(event1$customers, big.mark = ',')," customers."),
      group = "June 7-9 (blue)"
    ) %>%
    addPolygons(
      data = st_zm(event2),
      color = "#252525",
      weight = 2,
      smoothFactor = 0.1,
      fillOpacity = 0.7,
      fillColor = "#4daf4a",
      stroke = TRUE,
      label = ~event2$event,
      popup = paste0("This PSPS event spanned ", round(event2$Duration_hrs,0), " hours   \n   (beginning ", as.character(event2$start_time), " and ending at ", event2$end_time,").     \n   Ultimately PG&E estimates that it affected ",format(event2$customers, big.mark = ',')," customers."),
      group = "Sep 25-27 (green)"
    ) %>%
   addPolygons(
      data = st_zm(event3),
      color = "#252525",
      weight = 2,
      smoothFactor = 0.1,
      fillOpacity = 0.7,
      fillColor = "#984ea3",
      stroke = TRUE,
      label = ~event3$event,
      popup = paste0("This PSPS event spanned ", round(event3$Duration_hrs,0), " hours   \n   (beginning ", as.character(event3$start_time), " and ending at ", event3$end_time,").     \n   Ultimately PG&E estimates that it affected ",format(event3$customers, big.mark = ',')," customers."),
      group = "Oct 5 (purple)"
    ) %>%
   addPolygons(
      data = st_zm(event4),
      color = "#252525",
      weight = 2,
      smoothFactor = 0.1,
      fillOpacity = 0.7,
      fillColor = "#ffff33",
      stroke = TRUE,
      label = ~event4$event,
      popup = paste0("This PSPS event spanned ", round(event4$Duration_hrs,0), " hours   \n   (beginning ", as.character(event4$start_time), " and ending at ", event4$end_time,").     \n   Ultimately PG&E estimates that it affected ",format(event4$customers, big.mark = ',')," customers."),
      group = "Oct 9-12 (yellow)"
    ) %>%
       addPolygons(
      data = st_zm(event5),
      color = "#252525",
      weight = 2,
      smoothFactor = 0.1,
      fillOpacity = 0.7,
      fillColor = "#a65628",
      stroke = TRUE,
      label = ~event5$event,
      popup = paste0("This PSPS event spanned ", round(event5$Duration_hrs,0), " hours   \n   (beginning ", as.character(event5$start_time), " and ending at ", event5$end_time,").     \n   Ultimately PG&E estimates that it affected ",format(event5$customers, big.mark = ',')," customers."),
      group = "Oct 23-25 (brown)"
    ) %>%
    addPolygons(
      data = st_zm(event6),
      color = "#252525",
      weight = 2,
      smoothFactor = 0.1,
      fillOpacity = 0.7,
      fillColor = "#e41a1c",
      stroke = TRUE,
      label = ~event6$event,
      popup = paste0("This PSPS event spanned ", round(event6$Duration_hrs,0), " hours   \n   (beginning ", as.character(event6$start_time), " and ending at ", event6$end_time,").     \n   Ultimately PG&E estimates that it affected ",format(event6$customers, big.mark = ',')," customers."),
      group = "Oct 26 - Nov 1 (red)"
    ) %>%
       addPolygons(
      data = st_zm(event7),
      color = "#252525",
      weight = 2,
      smoothFactor = 0.1,
      fillOpacity = 0.7,
      fillColor = "#ff7f00",
      stroke = TRUE,
      label = ~event7$event,
      popup = paste0("This PSPS event spanned ", round(event7$Duration_hrs,0), " hours   \n   (beginning ", as.character(event7$start_time), " and ending at ", event7$end_time,").     \n   Ultimately PG&E estimates that it affected ",format(event7$customers, big.mark = ',')," customers."),
      group = "Nov 20-21 (orange)"
    ) %>%
    addLayersControl(
    # baseGroups = c("Street Map", "Terrain", "Sattelite", "B/W"),
      overlayGroups = c("HPI","June 7-9 (blue)", "Sep 25-27 (green)", "Oct 5 (purple)","Oct 9-12 (yellow)", "Oct 23-25 (brown)" , "Oct 26 - Nov 1 (red)",  "Nov 20-21 (orange)"),
      options = layersControlOptions(collapsed = FALSE)
                )%>%  # adds a button to return to the whole state view
      addLegend("bottomleft",
                pal = pal2,
                values = na.omit(maphpi$hpi2_pctile), 
                labels = c("bottom 25%","","","healthiest 25%"),
                opacity = 1,
                title = "HPI Percentile") %>%
      setView(lat = 39.302394, lng = -122.084003, zoom = 6) %>% #defaults the view of the original map to show the entire state
      addEasyButton(easyButton(
        icon="fa-globe", title="Zoom to State",
        onClick=JS("function(btn, map){ map.setView([37.085206, -119.540085],5); }"))) %>%
      hideGroup(c("June 7-9 (blue)", "Sep 25-27 (green)", "Oct 5 (purple)","Oct 9-12 (yellow)", "Oct 23-25 (brown)" , "Oct 26 - Nov 1 (red)",  "Nov 20-21 (orange)"))
    # adds a button to return to the whole state view
      # addLegend("bottomleft",
      #           pal = pal,
      #           values = factor(c("Moderate","High","VeryHigh"), levels = c("Moderate","High","VeryHigh")),
      #           opacity = 1,
      #           title = "Fire Hazard Severity Zones"
      #           ##lables = c("X%-X% (most vulnerable)","X%-X%","X%-X%","X%-X%","X%-X% (least"%. In this tract, the "mapTemp$def," is higher than "mapTemp),))##
      # )

Using the Health Equity Index

The census tracts impacted by PSPS can be quickly analyzed using existing Health Equity data from the Healthy Places Index (HPI).

The mean HPI of locations impacted by PSPS tend to be higher (healthier) than those not impacted.

foo <- merge({
  master_PSPS_tract[, .(shutoff = ifelse(PSPS_yn == "Yes", 1,0)), by = (GEOID) ] %>% .[, .(PSPS_yn = max(shutoff)), by = (GEOID) ] %>% .[, .(geoid = GEOID, PSPS_yn = ifelse(PSPS_yn ==0,"Not Affected by PSPS", "Affected by PSPS"))]
}, hpi) 

foo <- foo[,.(ct10, hpi2score, economic, education, housing, healthcareaccess, neighborhood, pollution, transportation, social, PSPS_yn)] %>% 
  melt.data.table(id.vars = c("ct10","PSPS_yn"), variable.name = "hpi_component",value.name = "value")


  hcboxplot(var = foo$hpi_component, x = foo$value, var2 = foo$PSPS_yn,
          outliers = FALSE) %>% 
  hc_chart(type = "column") %>% 
  hc_add_theme(hc_theme_smpl()) %>% 
  hc_xAxis(title = list(text = "HPI Index and Components")) %>%
  hc_yAxis(title = list(text = "HPI Score (higher = healthier)")) %>%
  hc_title(text = "Healthy Places Index and PSPS-affected areas",
           margin = 20, align = "left") %>%
  hc_subtitle(text = "data from PG&E and Public Health Alliance of Southern CA",
              align = "left", style = list(fontWeight = "bold"))

When we examine the number of times a place has been impacted, it appears that the places impacted the most in 2019 are worse off than those that were impacted once or twice, and in some cases worse than those places not impacted by PSPS.

This result is in part a function of the size of the area and population affected by the different PSPS events (larger events captured more areas that were only impacted by 1 or 2 PSPS events). We would like to examine the idea presented here using duration (hrs) of impact rather than frequency (#), particularly since different events last for different amounts of time. To do this we will need to improve the precision and accuracy of the data from PG&E.

foo <- merge({
  master_PSPS_tract[, .(shutoff = ifelse(PSPS_yn == "Yes", 1,0)), by = (GEOID) ] %>% 
    .[, .(PSPS_yn = max(shutoff), PSPS_count = sum(shutoff)), by = (GEOID) ] %>% 
    .[, .(geoid = GEOID, PSPS_count = ifelse(PSPS_count >=5, "5 or more",PSPS_count), PSPS_yn = ifelse(PSPS_yn ==1, "Affected by PSPS","Not Affected by PSPS"))]
}, hpi) 

foo <- foo[,.(ct10, hpi2score, economic, education, housing, healthcareaccess, neighborhood, pollution, transportation, social, PSPS_yn, PSPS_count)] %>% 
  melt.data.table(id.vars = c("ct10","PSPS_yn","PSPS_count"), variable.name = "hpi_component",value.name = "value")


  hcboxplot(var = foo$hpi_component, x = foo$value, var2 = foo$PSPS_count,
          outliers = FALSE) %>% 
  hc_chart(type = "column") %>% 
  hc_add_theme(hc_theme_smpl()) %>% 
  hc_xAxis(title = list(text = "HPI Index and Components")) %>%
  hc_yAxis(title = list(text = "HPI Score (higher = healthier)")) %>%
  hc_title(text = "Healthy Places Index and PSPS-affected areas",
           margin = 20, align = "left") %>%
  hc_subtitle(text = "data from PG&E and Public Health Alliance of Southern CA",
              align = "left", style = list(fontWeight = "bold"))

We could take a similar look at the data using each of the indicators in the Climate Change and Health Vulnerability Indicators Dataset (note that for CCHVIs high values indicate more vulnerability or worse off conditions).

foo <- merge({
  master_PSPS_tract[, .(shutoff = ifelse(PSPS_yn == "Yes", 1,0)), by = (GEOID) ] %>% 
    .[, .(PSPS_yn = max(shutoff)), by = (GEOID) ] %>% 
    .[, .(COUNTYFI_1 = GEOID, PSPS_yn = ifelse(PSPS_yn ==0,"Not Affected by PSPS", "Affected by PSPS"))]
}, cchvi_tracts[race == "Total" & latest =="Y"], by = "COUNTYFI_1") 



foo <- foo[,.(ct10, est, ind_strt, PSPS_yn)] #%>% 
  #melt.data.table(id.vars = c("ct10","PSPS_yn"), variable.name = "hpi_component",value.name = "value")


  hcboxplot(var = foo$ind_strt, x = foo$est, var2 = foo$PSPS_yn,
          outliers = FALSE) %>% 
  hc_chart(type = "column") %>% 
  hc_add_theme(hc_theme_smpl()) %>% 
  hc_xAxis(title = list(text = "CCHIV Indicators")) %>%
  hc_yAxis(title = list(text = "CCHVI Values (higher = more vulnerable)")) %>%
  hc_title(text = "Climate Change & Health Vulnerability Indicators and PSPS-affected areas",
           margin = 20, align = "left") %>%
  hc_subtitle(text = "data from PG&E and CDPH",
              align = "left", style = list(fontWeight = "bold"))

Statistical Testing of Specific Indicators

Each of the plots above are boxplots, which describe the distribution of each indicator in the different areas. The center line in those plots represent the median for the group, and the whiskers are extremes. Moving forward, specific indicators (by event) could be investigated in a more rigorous way using data from the American Community Survey (ACS) directly to compare the estimated indicator values for different areas and to provide confidence intervals, which allow for statistical testing of differences between PSPS-affected areas and non-affected areas.

Below is an example for the % of people over 65 in the different areas by event. Similar comparisons have also been done for the number of adults living with a disability (and below the poverty line).

old <- tidycensus::get_acs(geography = "tract", 
              variables = c(numerator = "B01001_020", 
                            numerator = "B01001_021",
                            numerator = "B01001_022",
                            numerator = "B01001_023",
                            numerator = "B01001_024",
                            numerator = "B01001_025",
                            numerator = "B01001_044", 
                            numerator = "B01001_045",
                            numerator = "B01001_046",
                            numerator = "B01001_047",
                            numerator = "B01001_048",
                            numerator = "B01001_049",
                            denominator = "B01001_001"), 
              state = "CA") %>% data.table()

old <- old[,.(estimate = sum(estimate), 
              moe = tidycensus::moe_sum(estimate = estimate, moe = moe)), 
           by= .(GEOID, NAME, variable)]

try1 <- merge(old, {
  readRDS("//mnt/projects/ohe/PSPS/PSPS_yn_master_tract.RDS")
  }, allow.cartesian = T)
# try1 <- na.omit(try1)

try1 <- try1[, se:= moe/1.645]

sums <- try1[,.(sum = sum(estimate, na.rm=T), 
                se_sum = (sum(se^2, na.rm=T)^0.5)), by = .(variable, PSPS_yn, event)]

totals <- sums %>% dcast.data.table(formula = PSPS_yn + event ~ variable, value.var = "sum")

sum_ses <- sums %>% dcast.data.table(formula = PSPS_yn + event ~ variable, value.var = "se_sum") %>% .[,.(PSPS_yn, event, SE_num = numerator, SE_pop = denominator)]

final <- merge(totals, sum_ses)

final <- final[, `:=` (proportion = numerator/denominator,
                       se_prop = (1/denominator)*(SE_num^2 - (((numerator/denominator)^2)*SE_pop^2))^0.5)]

final <- final[, `:=` (LL = proportion - 1.96*se_prop,
                       UL = proportion + 1.96*se_prop)]

final <- final[, event := factor(event, levels = c("June_7-9", "September_25-27","October_5","October_9-12","October_23-25","October_26-November_1","November_20-21"))] %>% setorder(event)

final <- final[, PSPS_yn := ifelse(PSPS_yn == "No","Not Affected by PSPS", "Affected by PSPS")]



hc1 <- highchart() %>% 
  hc_chart(type = "bar") %>%  
  hc_add_series(final, "point", hcaes(y = 100*proportion, x = event, group = PSPS_yn))%>% 
  hc_add_series(final, "errorbar", hcaes(low = 100*LL, high = 100*UL, x = event, group = PSPS_yn)) %>%
  hc_add_theme(hc_theme_smpl()) %>% 
  hc_xAxis(title = list(text = "PSPS Event"), 
           type = 'category'  ,
          categories =  c("June 7-9  (14hrs, 22k)", "Sep 25-27  (26hrs, 75k)","Oct 5  (18hrs, 11k) ","Oct 9-12  (90hrs, 729k)","Oct 23-25  (52hrs, 177k)","Oct 26-Nov 1  (129hrs, 941k)","Nov 20-21  (42hrs, 49k)")
          ) %>%
  hc_yAxis(title = list(text = "% over 65 (with 95% CI)")) %>%
  hc_title(text = "% of population over 65",
           margin = 20, align = "left") %>%
  hc_subtitle(text = "Tracts Impacted by PSPS vs Not Impacted",
              align = "left", style = list(fontWeight = "bold"))


hc1

This Table summarizes much of the data we currently have from PG&E on their events as well a few of the comparisons on indicators of vulnerability for areas impacted and not impacted in each event.

Jun 7-9 Sep 25-27 Oct 5 Oct 9-12 Oct 23-25 Oct 26-Nov 1 Nov 20-21
start time 6/8/2019 6:00 9/23/2019 17:06 10/5/2019 22:00 10/9/2019 0:09 10/23/2019 13:54 10/26/2019 16:00 11/20/2019 3:38
end time 6/8/2019 20:00 9/24/2019 18:40 10/6/2019 16:00 10/12/2019 17:41 10/25/2019 18:20 11/1/2019 1:20 11/21/2019 21:55
Duration_hrs 14 25.57 18 89.53 52.43 129.33 42.28
customers 22,327 75,385 11,304 728,980 176,620 941,217 49,085
residential 19,500 35,501 9,981 636,355 832,314
comm_indust 2,565 3,856 1,250 81,318 98,514
med_baseline 1,589 4,451 718 30,026 7,823 34,618 2,456
other 262 616 73 11,307 10,389
% Over 65
Affected by PSPS 22.8% (22.1, 23.5) 23.6% (22.8, 24.4) 24.3% (23, 25.7) 17.9% (17.7, 18) 24.5% (24, 25) 18.7% (18.6, 18.9) 24.2% (23.3, 25.1)
Not Affected 13.2% (13.1, 13.2) 13.2% (13.1, 13.2) 13.2% (13.2, 13.2) 12.9% (12.8, 12.9) 13.1% (13, 13.1) 12.7% (12.7, 12.8) 13.2% (13.1, 13.2)
% Disabled
Affected by PSPS 11.1% (10.5, 11.8) 19.6% (18.7, 20.5) 23.2% (21.5, 24.9) 11.9% (11.8, 12.1) 15.9% (15.4, 16.4) 11.6% (11.5, 11.8) 16.9% (15.8, 17.9)
Not Affected 9.9% (9.9, 10) 9.9% (9.8, 9.9) 9.9% (9.9, 9.9) 9.8% (9.7, 9.8) 9.8% (9.8, 9.9) 9.8% (9.7, 9.8) 9.9% (9.9, 9.9)
% Disabled and Below Poverty Line
Affected by PSPS 1.5% (1.3, 1.8) 3.8% (3.4, 4.2) 4.4% (3.6, 5.2) 2.1% (2, 2.2) 2.5% (2.3, 2.8) 1.9% (1.8, 2) 3.2% (2.7, 3.7)
Not Affected 2% (2, 2) 2% (1.9, 2) 2% (2, 2) 2% (1.9, 2) 2% (1.9, 2) 2% (2, 2) 2% (2, 2)